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ABSTRACT 

We demonstrate the existence of an enhanced rate of angular momentum 
relaxation in nearly Keplerian star clusters, such as those found in the centers of 
galactic nuclei containing massive black holes. The enhanced relaxation arises 
because the radial and azimuthal orbital frequencies in a Keplerian potential are 
equal, and hence may be termed resonant relaxation. We explore the dynamics 
of resonant relaxation using both numerical simulations and order-of-magnitude 
analytic calculations. We find that the resonant angular momentum relaxation 
time is shorter than the non-resonant relaxation time by of order M*/M, where 
M* is the mass in stars and M is the mass of the central object. Resonance 
does not enhance the energy relaxation rate. We examine the effect of resonant 
relaxation on the rate of tidal disruption of stars by the central mass; we find 
that the flux of stars into the loss cone is enhanced when the loss cone is empty, 
but that the disruption rate averaged over the entire cluster is not strongly 
affected. We show that relativistic precession can disable resonant relaxation 
near the main-sequence loss cone for black hole masses comparable to those 
in galactic nuclei. Resonant dynamical friction leads to growth or decay of 
the eccentricity of the orbit of a massive body, depending on whether the 
distribution function of the stars is predominantly radial or tangential. The 
accelerated relaxation implies that there are regions in nuclear star clusters 
that are relaxed in angular momentum but not in energy; unfortunately, these 
regions are not well- resolved in nearby galaxies by the Hubble Space Telescope. 



Subject headings: black hole physics — galaxies: active — galaxies: kinematics 
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The dynamical evolution of galactic and active galactic nuclei (AGNs) is determined 
by numerous physical processes, including two-body relaxation, stellar collisions/mergers, 
tidal disruption, binary formation, stellar evolution, and disk hydrodynamics (e.g., Spitzer 
& Stone 1967; Norman & Scoville 1988; Quinlan & Shapiro 1990; Murphy et al. 1991). In 
most cases, however — and in contrast to systems such as globular clusters — the two-body 
relaxation times are so long that relaxation has had little observable effect on the nucleus 
so far. The slow relaxation rates in galactic nuclei are responsible, for example, for the 
general failure of stellar tidal disruption by relaxation into the loss cone (Frank & Rees 
1976; Lightman & Shapiro 1977) as a viable AGN fueling mechanism (Hills 1975; Frank 
1978; McMillan et al. 1981), and for the absence of core collapse or mass segregation at 
observable resolutions in all but a few nearby galaxies (e.g., M33; Kormendy & McClure 
1993). 

Normally the gravitational potential in a galaxy is determined by the stars and 
distributed dark matter. However, in galactic nuclei possessing massive black holes, the 
potential in the inner nucleus is dominated by the black hole and hence nearly Keplerian, 
so that eccentric orbits maintain their spatial orientation for many orbital periods. The aim 
of this paper is to demonstrate that in Keplerian and other such "resonant" potentials the 
rate of relaxation of angular momentum can be greatly enhanced, so that relaxation can be 
important even if the energy relaxation timescale is longer than a Hubble time. We examine 
the dynamics of this "resonant relaxation" through approximate analytic arguments and 
iV-body simulations, and discuss its influence on galactic nuclei. 

1. The Process of Resonant Relaxation 
1.1. Introduction 

The force field F(r, t) in an equilibrium iV-body stellar system can be divided into a 
mean force F(r) = (F(r,t)} and a fluctuating force f(r,t) = F(r,t) — F(r), where (•) 
denotes time average. If N 3> 1 the fluctuating force is a Gaussian random field, and hence 
is completely described by the correlation function Cij(ri,r2,r) = (fi(vi,t), fj(r2,t + r)). 
This fluctuating force induces diffusion or relaxation of the deterministic orbits that stars 
would follow if only the mean force field F were present. The relaxation time, t re \, is 
crudely defined so that the diffusion of integrals of motion such as energy E and angular 
momentum L (per unit mass) is given by AE ~ E (t/trci) 1 ^ 2 and AL ~ L (t/t^i) 1 ^ 2 . 

The usual estimate of the relaxation time (Jeans 1913, 1916; Chandrasekhar 1942; 
Binney & Tremaine 1987) is based on an infinite homogeneous stellar system in which stars 
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travel on straight-line orbits. This assumption is plausible because each octave in spatial 
scale between the system size R and the much smaller scale r min ~ R/N (the scale on which 
close encounters produce ~ 90° deflections) contributes equally to the relaxation rate, and 
for most of these octaves the approximations of homogeneity and straight-line orbits are 
legitimate (summing the contributions from different scales gives rise to the well-known 
Coulomb logarithm In A ~ ln(i?/r min ) ~ ln(iV) which appears in formulae for the relaxation 
rate). One consequence of the assumption of straight-line orbits is that the correlation 
function Cy(ri, r 2 , r) decays rapidly to zero when r exceeds \r\ — r 2 \/V, where V is a 
typical velocity (for infinite homogeneous systems — > 1/r as r — > oo [e.g. Cohen 1975]). 
Our focus here is on stellar systems in which the correlation function remains non-zero for 
much larger times, a condition that occurs if most of the stars are near resonance. 

If motion in the mean force field F(r) is regular, then stellar orbits are quasiperiodic 
with three characteristic frequencies f^. The orbits are resonant if there are linear 
combinations of the form X^=i = where the h L are small integers. The simplest 
important examples are (i) spherical potentials, in which one frequency is zero because all 
orbits remain in a fixed plane; (ii) Kepler potentials, in which one additional frequency is 
zero because the apsis does not precess; and (iii) the harmonic oscillator potential, in which 
the radial frequency is twice the azimuthal frequency, so that the orbit shape is a centered 
ellipse. 

The possibility of enhanced relaxation in potentials that support many near-resonant 
orbits was discussed by J. Ostriker two decades ago, in lectures for a graduate course in 
stellar dynamics attended by one of us (Ostriker 1973). 



1.2. Non-resonant Relaxation 

We begin by examining relaxation in a near-Kepler potential. Consider a spherical 
volume of radius R centered on a point mass M and containing iV ^> 1 identical stars 
of mass m, where M* = Nm <C M. We assume that the stellar orbits have random 
orientations and moderate eccentricities, and that the density of stars is approximately 
uniform within R. The typical stellar velocity is V ~ (GM/R) 1 / 2 and the characteristic 
orbital period is t orb ~ R/V. Since <C M, each orbit is approximately a Kepler ellipse, 
which precesses slowly on a timescale t pre c- If the precession is dominated by the mean field 
from the other stars (rather than, say, relativistic effects or an external tidal field), then 

M 

tprec ~ ~^~torhi (1) 

which is much longer than t or b. 
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The usual (non-resonant) relaxation rate can be estimated by the following argument. 
Consider a volume of radius r < R, which typically contains n ~ N (r/R) 3 stars. The 
instantaneous number of stars in this volume fluctuates by an amount n 1//2 ; thus the 
fluctuating force on the length scale r is / ~ Gmn 1 ^ 2 jr 2 . The force fluctuates on a 
timescale tf ~ r/V, and the typical impulse that a star receives during this timescale 
is Sv ~ ftf ~ Gmn 1 / 2 / (rV). Since impulses in successive intervals of length tf are 
uncorrelated, the total impulse after time t is described by a random walk, 

. . , 9 t G 2 m 2 n G 2 m 2 N 9 m 2 N t 

Note that At> is independent of the scale r, which confirms that each octave in scale 
contributes equally to the relaxation rate; thus the total diffusion rate must be multiplied 
by a factor In A which represents the number of octaves that contribute to the relaxation. 
We may define the non-resonant relaxation time by (Av/V) 2 = (t/t™\), so that 

M 2 

^rel ~ 2~7\M A~ ^orb- (3) 

m z i\ in A 

The fluctuating force changes both energy and angular momentum, at rates given by 

mN 1 / 2 

(AE/E) nT ~ (*/C) 1/2 ~ «^-(iAo rb ) 1/2 , 

771 AT 1 / 2 

(AL/L max ) nr ~ (i/Ci) 1/2 ~%^-(iAorb) 1/2 , (4) 

where E ~ V 2 ~ GM/R, i^ax ~ ^ 2 -R 2 ~ GMR, and ct and r] s are dimensionless constants 
that equal the square root of the Coulomb logarithm to within a factor of order unity. 



1.3. Resonant Relaxation in Near-Kepler Potentials 

To estimate the resonant relaxation rate, we imagine averaging the stellar density over 
an intermediate timescale that is ^> t OTh but t prC c- On this timescale each star can be 
represented by a fixed wire whose mass is the stellar mass, whose shape is a Kepler ellipse, 
and whose linear density is inversely proportional to the local speed in the elliptical orbit. 

The gravitational potential from these wires is stationary and hence does not lead to 
energy relaxation; thus 

(AE/E) ICS = 0. (5) 

However, the wires exert mutual torques which induce angular momentum relaxation. The 
typical specific torque on a wire is T ~ N l l 2 Gm/ R, and fluctuates on a timescale ~ t prcc as 
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the wires precess in different directions. Thus the characteristic change in specific angular 
momentum AL ~ Tt over a timescale t < t pTec is given by 

mN 1/2 

( A L I Z/ ma x ) res ~ Ps M (t/t OTh ), t<tprec, (6) 

where f3 s is a dimensionless constant of order unity. 

Over timescales t 3> t pr ec, the change in angular momentum is described by a random 
walk with increments AL/L max ~ ^(miV 1 / 2 /M)(t pTec /t OJ .y ) ) over a characteristic time £ pre c; 
thus, 

fKTlr \ a mNl/2 fWA 1/J f>7\ 

(AL/L max ) rcs ~ p g — — — \ —p — I , t » t prec ; (7) 

if the precession is determined by the mean field of the other stars then equation ([I]) implies 
that 



(Tfl \ 
Tt) (^Aorb) 1//2 , ^ > Ve - 



in this case the angular momentum relaxation time, defined by (AL/L max ) res ~ (t/t^f) 1 / 2 , 
is given by 

M 

^rel ~ *orb, (9) 

m 

which, remarkably, is independent of the number of stars N. The resonant relaxation 
time is shorter than the non-resonant relaxation time (eq. 0) by a factor (mN/M) In A. 
If there is a range of stellar masses, the factor m in equation (§) should be replaced by 
/ m 2 dN(m)/ f mdiV(ra). 

There is no analog to the Coulomb logarithm in resonant relaxation, since the 
relaxation is dominated by large-scale fluctuations. The absence of the Coulomb logarithm 
implies that the resonant diffusion rate depends on the overall structure of the stellar 
system and cannot be computed using the assumption of local homogeneity, as is done for 
the non-resonant relaxation rate. 

On timescales longer than the resonant relaxation time but shorter than the 
non-resonant relaxation time, a stellar system in a Kepler potential should be in the 
maximum-entropy state consistent with its original total angular momentum L tot and 
energy distribution N(E) (the number of stars with energy < E, which is invariant on 
this timescale since there is no resonant energy relaxation). This state is described by the 
phase-space distribution function 



f(r,v) = w(E) exp(— b • L), 



(10) 



- 6- 



1„,2 
2 

constraints 



where E = ^v 2 — GM/r, L = r x v, and b and are determined implicitly by the 

/ f(r, v)L drdv, 

(11) 

J f(r 1 ,v 1 )S(E-E 1 )dr 1 dv 1 . 



Ltot 

dN(E) 
dE 



If L to t = then 6 = and the distribution function is a function of energy alone (i.e., the 
cluster is isotropic). If there is a range of masses then (flO|) is replaced by 

f(r,v,m)=w(E,m)exp(—mb-L). (12) 



1.4. Resonant Relaxation in Near-spherical Potentials 

A more limited form of resonant relaxation is present in any spherical potential. 
Consider a generic spherical potential and average the stellar density over a timescale 
^> t or b. On this timescale each star is smeared into an axisymmetric annulus whose inner 
and outer radii are the pericenter and apocenter distances. The gravitational potential 
from these annuli is stationary and hence does not lead to energy relaxation; thus, as in 
equation (|5|), 

{AE/E) ies = 0. (13) 

The annuli exert mutual torques which induce angular momentum relaxation; however, in 
contrast to the Kepler case the torques are perpendicular to the orbit normals (because 
the averaged orbits are axisymmetric annuli rather than eccentric wires). Thus the 
vector torque between two annuli with vector angular momenta Li and Lj satisfies 
Tij ■ Li = Tij ■ Lj = 0; in other words the torques change the directions of the angular 
momentum vectors but not their magnitudes. Thus, in contrast to equation (|]), there is no 
resonant relaxation of the scalar angular momentum, 

(AL/L max ) rcs = (A|L|/L max ) res = 0, (14) 

but there is resonant relaxation of the vector angular momentum. The resonant relaxation 
rate may be estimated by analogy with equation (||). The typical specific torque on 
an annulus is T ~ N x l 2 Gmj R, and fluctuates on a timescale ~ £p rec - Here t^ mc is the 
precession time for the angular momentum vector (not the apsis, as in §§ |1.2j and 
is determined by the stochastic component of the potential, t^ rec ~ A^ 1 / 2 t or b//i. 

The characteristic change in vector angular momentum over a timescale t < t^ rec is 
given by \AL\ ~ Tt. Since we are considering general near-spherical potentials, we no 
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longer require that the potential is dominated by a black hole, M ± <C M; thus the maximum 
angular momentum is given by 1^ ~ V 2 R 2 ~ G(M + M*)R and we have 



(|A£|/L max ) rcs ~ /i A^(t/t orb ), t < t^ rec , (15) 



where /3 V is a dimensionless constant of order unity, and 

^a^Tm (16) 

(~ 1 if the potential is dominated by the stars themselves). 

Over timescales t 3> t^ Tec , the change in angular momentum is described by a random 
walk. The increments AL of the random walk are determined by evaluating (|15|) at t ~ £p rec , 
which yields \AL\/L maiX ~ 1. In other words the angular momentum vectors drift over the 
whole velocity sphere on a timescale of order tp rec , implying t^i ~ t^ vec ~ A^ 1 / 2 t orb / / u, which 
is shorter than the non-resonant relaxation time t™ { (eq. 0) by a factor (fi/N 1 ^ 2 ) InA. 

A closely related form of resonant relaxation is present in general axisymmetric 
potentials that are nearly spherical. Once again there is no resonant relaxation of the scalar 
angular momentum, but there is resonant relaxation of the vector angular momentum. In 
this case only the z-component of the vector angular momentum is conserved by motion in 
the mean potential, so we focus on this component. By analogy with equation ( |15|) we have 

(|AL 2 |/L max ) res ~ ^-^(t/t orb ), t « ^ rec , (17) 

where in this case tp rec , the precession time for the angular momentum vector, is determined 
by the non-spherical component of the mean potential. Over timescales 2> £p rec , the change 
in L z is described by a random walk with increments |AL 2 |/L max ~ (A i /^z/^ 1 / 2 )(tp I . ec /t rb) 
over a characteristic time £p rec ; thus, 

(|AL 2 |/L max ) rcs ~ ^ 7 , t » t L pmc . (18) 

The relaxation time, defined by (| AL 2 |/L max ) res ~ (^Arei) ; i s then given by 

N t 2 u 

fc rel 2 fL ' 

h 1 b prcc 

which is shorter than by a factor (t rb/^p re c) ^ n -^- 

The vector angular momentum also diffuses from non-resonant relaxation, at a rate (cf. 
eq. 1) 

(|AX|/L max ) nr ~ /i^(t/tcrb) 1/a , (20) 
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where r] v is a dimensionless constant that equals the square root of the Coulomb logarithm 
to within a factor of order unity. 

1.5. Resonant Relaxation in Regular Potentials 

Motion in any time-independent, regular potential can be described by action-angle 
variables (J*, fa) and a Hamiltonian H(J). For spherical potentials the actions (Ji, J 2 , J3) 
can be chosen to be, respectively, the radial action, the total angular momentum, and the 
^-component of the angular momentum (e.g., Tremaine & Weinberg 1984). The motion is 
quasiperiodic with fundamental frequencies fa = = dH/dJi. In a spherical potential, 

Qi is the radial frequency, fl 2 is the azimuthal frequency, and Q3 = 0; in a Kepler potential 
Qi — Q'2- Resonant relaxation can be important if the resonance condition Yli=i k&i = 
is approximately satisfied for most stars, where the fcj's are small integers with no common 
factor (e.g., k\ = k 2 = for a spherical potential; k\ = —k 2 for a Kepler potential). 

To isolate the effects of resonant relaxation we perform a canonical transformation to 
new action-angle variables (Ki,ipi) defined by the generating function 

3 

S(K h fa) =KxY, hfa + K 2 fa + K 3 fa. (21) 

i=l 

Then 

Jx = k x K x , J 2 = k 2 K 1 + K 2 , J 3 = k 3 K 1 + K 3 , (22) 

and 

3 

^1=^2 kifa, ^ 2 = 02, ^3 = 03- (23) 

1=1 

The resonance condition implies that ipi is slowly varying; the effects of resonant relaxation 
are therefore described by a fluctuating Hamiltonian of the form h(K,ipi). The resonant 
Hamiltonian does not depend on the non-resonant angles ip 2 and ipz or on the time, so 
Hamilton's equations imply that the conjugate momenta K 2 and K 3 and the total energy 
E = H + h are unaffected by resonant relaxation; this in turn implies that the changes in 
the energy and actions caused by resonant relaxation satisfy the constraints 

AE = 0, AJ = C(t)k, (24) 

where C(t) is scalar. 

Resonant relaxation leads to diffusion of C(t) which is described by analogs of equations 
(0) and 

(AC/L max ) res ~ ft , j-y 12 (^Aorb) 5 t <C ^p rec i 
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1/2 (25) 

(AC/L max ) rcs ~ ( "pT^ - J ' f » tc 

where 7 is a dimensionless constant of order unity, // is defined in equation fllED, 

tp rcc ~ (J2 h^i)' 1 is the time required for the slow angle to change by one radian, and 

we have assumed that [ As I is small. 



1.6. Resonant Friction 

Chandrasekhar (1943) showed that the stochastic changes in velocity associated with 
relaxation are accompanied by a systematic drag which he named "dynamical friction." 
For massive objects (MOs) orbiting in a stellar system, orbital evolution from dynamical 
friction is much faster than evolution from stochastic relaxation. Near-resonances can 
enhance dynamical friction just as they enhance stochastic relaxation, leading to an effect 
we may call "resonant friction." 

To analyze resonant friction we use the formalism developed by Lynden-Bell & Kalnajs 
(1972); the z-component of the specific torque on the MO from dynamical friction is 
(Tremaine & Weinberg 1984, eq. 65) 

T z =4^m ]T k 3 f dJ^ki^jl^kfSik-n-cUk). (26) 

fc,fc 3 >0 J i aJi 

Here mo > m is the mass of the MO, /(«/) is the phase-space density of bound stars, and 
the potential of the MO has been written in the form 

U(r,t) = m Re| ]T ^k(J) exp[i(k ■ 4> - to k t)} \ . (27) 



Equation ( ^6|) contains both resonant and non-resonant contributions to dynamical friction. 



Now let us assume that the stellar system is spherically symmetric and that the 
distribution function depends only on energy and angular momentum, / = f(E,L), 
consistent with Jeans' theorem. Without loss of generality we may assume that z = is 
the orbital plane of the MO, so that the total specific torque on the MO is T = T z . Using 
OE/dJi = dH/dJi = fij and L = J 2 we have 

T = 4vr 4 m £ k 3 [ dJ L k % + kM W fc | 2 5(k ■ tl - u k ). (28) 
k,k 3 >o J V dh dL J 
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Now suppose that there is a near-resonance for the triplet k r . Formally, we set 

k r ■ Q = e k r ■ O and Uk r = euJk r and let e — > 0. Keeping only the largest terms in fl2"8|) gives 

the contribution of the triplet k r to the specific torque from resonant friction, 

which diverges as e —>■ 0. 

Evaluating expressions such as (^9|) is arduous, mostly because of the complexity of 
the action-angle expansion of the potential of a point mass (Weinberg 1986; see also 
Hernquist & Weinberg 1989). However, it is simple to derive some qualitative effects of 
resonant friction: 



Equation ( |29| ) shows that there is no resonant friction when the distribution function 
is isotropic (df/dL = 0). The term oc df/dE in equation (p8|) contributes only 
non-resonant friction. 

The resonant torque is zero if the stars are on circular orbits (in this case = 
unless k\ = 0, and if ki = then resonance requires k 2 = since f2 2 7^ and Q 3 = 0); 
similarly, the resonant torque is zero if the MO is on a circular orbit. 

The orbits that contribute most to the friction are those with near-zero inclinations, 
since these remain closest to the zero-inclination orbit of the MO and precess in 
the same direction. Zero-inclination orbits have \I/fc = unless k 2 = k% (Weinberg 
& Tremaine 1984). Together with equation (p9]), this suggests that the dominant 
contribution to the resonant torque comes from terms with k^^k2 )T > 0, which in turn 
implies that the sign of the torque on the MO is the sign of df/dL. In the common 
case where orbits are predominantly radial, df/dL < 0, resonant friction removes 
angular momentum from the MO orbit at constant energy, thereby increasing its 
eccentricity. This effect can dominate the eccentricity evolution of black-hole binaries 
in galactic nuclei and may promote the merger of binary black holes, since emission of 
gravitational radiation is more efficient for an eccentric binary (Begelman et al. 1980, 
Quinlan 1996). 

The rate of growth or decay of angular momentum of the MO through resonant 
friction is approximately given by 




Lmax \ dt J L m£LX Af 2 
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where we have assumed that the orbit of the MO has moderate eccentricity and that 
\df/dL\ ~ f / L max . In a near-Kepler potential, where k r = (1,-1,0) and t* rec is given 
by equation (JTJ), we have 

1 fiL ) (3D 



Lmax \ dt J res M t or b 

In this case the frictional torque is — remarkably — independent of the number of stars, 
although equation (|3T| ) holds only when m,n ^ M*. When mo ^ M* the precession 
time is t pTec ~ (M/mn)t rb, so that 

1 / dL\ M* 1 

~±^^- (32) 



£max \ dt y res ill t rb 

The two expressions (|3lD and (^) can be combined, 



1 ( dL\ ^ ± min (M*,m ) 1 



imax \ dt y res ikf t or b 

For comparison, non-resonant dynamical friction removes energy and angular 
momentum at the slower rate 

L(*E\ ~J-(**L\ „-I^JL (34 ) 
E \ dt ) m L max V dt J nr M 2 t olh V ; 

Another manifestation of resonant friction occurs when the potential is nearly spherical 
but the mean rotation of the stars is non-zero (for example a stellar disk) and the orbit 
of the MO is inclined. In this case the distribution function / = f(E, L, L z ) depends not 
just on E and L but also on the ^-component of angular momentum L z = J3. Because 
the potential is nearly spherical, Q3 ~ (i.e. inclined orbits precess slowly) so terms with 
k\ = k 2 = are near-resonant. The z-component of the torque on an MO from the triplets 
k = (0,0, k 3 ) is 

r df 2 

fc 3 >0 J 0Lz 

If the stellar system rotates in a prograde direction, then df/dL z is generally positive so 
T z is positive; thus the resonant friction erodes the inclination of the MO until it settles 
into the equatorial plane of the stellar system. The resonant torque is zero if the stars have 
precisely zero inclination (for zero- inclination stars ^/o,o,fc 3 = unless ks = 0); the rate of 
change of the orbital inclination In of the MO is given approximately by 



1 / dip \ M^mp t prcc 

In I dt J ~ M 2 t 2 



(I 2 ), (36) 



'0 V ut / res 1V1 b orb 

where (I 2 ) is the mean-square inclination in the disk and we have assumed 
\df/dL z \~f/L 

max • 
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2. Numerical Simulations 
2.1. Methods 

Our numerical investigation of resonant relaxation utilized two complementary 
approaches: an iV-body code that integrates the orbits of individual stars with timestep 
<C t or b, and an "iV-wire" code that follows the evolution of Kepler ellipses with timestep 
<C tprec- The iV-body code allows us to follow the growth of AE(t), AL(t), and AL(t) for 
each star; the iV-wire code involves a degree of abstraction from the original problem and 
does not yield AE(t), but potentially allows AL(t) and AL(t) to be followed for much 
longer times. 

N-body simulations These simulations are challenging because isolating the effects 
of resonant relaxation requires following the iV-body system for very long times (our 
simulations ran for up to 10 6 t or b). 

The iV-body system contained three components: 

• A fixed spherical potential, $(r). To measure the importance of resonant relaxation 
we compared the evolution of pairs of systems that were identical except for the fixed 
potentials. One potential was Keplerian, $(r) = —GM/r, in which there is no orbital 
precession, and the other was an isochrone potential, <£>(r) = —GM/[b + (b 2 + r 2 ) 1 / 2 ] 
(Henon 1959), where b is a "core" radius, in which the precession time is comparable 
to the orbital period for r ~ b. We used units in which G = M = b = 1. 

• N identical "background" stars of mass m, which felt only the gravitational force 
from the fixed potential, not from each other or from the test stars. The total mass 
of the background stars was much smaller than the mass associated with the fixed 
potential; typically M+/M ~ W^-IO' 5 . 

• n -C iV identical "test" stars of mass m, which felt the gravitational force from the 
central potential, the background stars, and from one another. (We also conducted a 
few simulations using test stars of unequal masses.) The use of separate background 
and test stars reduces the force calculation to 0(nN) rather than the much larger 
0(N 2 ) for the fully self-consistent calculation; typically, n/N ^ 0.1. 

To reduce integration errors during close encounters, the gravitational force between 
stars was softened. The softening length was usually ~ 1% of the system size; the influence 
of softening on the results is discussed in § |2.2j . 
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The initial distribution function of the stars was isotropic, with dN(E)/ dE oc E~ 2 
for E min < E < E max = 10-E min and zero otherwise. For the isochrone potential, 
l-E'minl = 0.05 GM/b; this implies a stellar density p(r) oc r~ 2 for b <^ r <^ 10b. For the 
Kepler potential, E min was chosen so that the typical star had the same orbital period as 
in the isochrone potential. Test stars were chosen such that their orbits remained within 
the body of the cluster, so as to avoid edge effects due to the cluster's finite radial extent; 
in practice this meant limiting possible test stars to orbits that were not too eccentric, 
e <^ 0.8. For example, in the isochrone potential test stars were required to have turning 
points r min > 0.8 6 and r max < 9 6, implying e < 0.84; background stars were present in the 
larger region 0.2 b < r < 20b. 

The integration algorithm was symplectic to eliminate spurious dissipation over 
the long integration times. We used either a second-order scheme ("leapfrog") or its 
fourth-order generalization^, depending on the situation — the higher-order scheme was 
more accurate but the lower-order scheme was more robust during close encounters. An 
additional advantage of these algorithms is that they exactly conserve angular momentum 
when integrating motion in a central potential. 

The advantages of symplectic schemes are generally lost if the timestep is varied. Thus 
each orbit was integrated with a fixed timestep. However, because of the wide range of 
orbital periods, it was useful to allow different timesteps for different orbits. To avoid having 
to interpolate when computing forces, the different timesteps must be commensurate; in the 
case of second-order leapfrog, this is easily arranged. Suppose that a timestep less than hi is 
needed to maintain the desired accuracy for star i, and let h = min^ hi. In a leapfrog scheme 
with step h the forces are evaluated at times (k+ ^)h, k = 0, 1, 2 . . ., that is, in the middle 
of the timestep. Therefore we can ensure that the timesteps are commensurate by choosing 
timesteps h[, where h[ is the largest odd multiple of h less than hi. The further restriction 
that h'i/h be a multiple of 3 ensures that all of the timesteps end at the same time as the 
largest timestep, making synchronized output easy. Using this approach speeded up the 
iV-body code by a factor ~ 2-3. In the fourth-order scheme, however, forces are computed 
at irrational (and hence unalignable) fractions of the interval, so that in this case the same 
timestep h was used for all stars; for fractional energy accuracies e^3 x 10 -5 , the savings 
of the fourth-order scheme outweighed the cost of a single timestep and produced the faster 
code. 



1 The leapfrog map with timestep h is H 2 {h) : (r,v) — > (r',v'), where r\ = r + ^hv 7 v 1 = v — /iV£/(ri), 
r' = r% + i/iu'; the fourth-order map (e.g., Yoshida 1990) is the composition of three leapfrog steps, 
H 4 {h) = H 2 (ah)H 2 {bh)H 2 (ah), where a = 1/(2 - 2 1 / 3 ) and b = 1 - 2a = -2 1 / 3 /(2 - 2 1 / 3 ). 
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The size of the timestep (in units of the radial orbital period) required to maintain 
relative energy accuracy e was h ~ (^min/^max) 1 ' 25 ^ 1 ^ 2 for the second-order scheme and 
h ~ 0.5(r min /r max ) 1 - 25 e 1 / 4 for the fourth-order scheme, where r min and r max are the orbit's 
turning points (e is the amplitude of the spurious energy oscillation caused by the integration 
algorithm; there is no systematic energy drift, because the algorithm is symplectic). A 
useful check on the accuracy of the code is to follow a cluster whose background stars are 
artificially held fixed at their initial positions; since the test stars are then orbiting in a fixed 
(albeit non-spherical) potential, their orbital energy should be conserved. This test was 
mainly used to ensure that the error accumulated during close encounters was negligible. A 
second constraint on the required accuracy is that the natural orbital precession should not 
be artificially accelerated by integration errors. The simulations generally used e = 10~ 6 . 

N-wire simulations The heart of these simulations is an efficient algorithm to compute 
the time-averaged torque between two Kepler orbits, which was supplied to us by Jihad 
Touma (Touma & Tremaine 1996). The evolving angular momentum for each test star was 
computed by straightforward integration of the time-dependent torque on the orbits using 
a standard (non-symplectic) adaptive ODE integrator. 

The initial conditions and the division into test and background stars were the same as 
in the A-body code. Of course, since the analysis assumes that the orbits are approximately 
Keplerian, comparative integrations in the isochrone potential were not possible. The 
simulations also do not yield estimates for AE, which is zero in this approximation. The 
potential advantages of the iV-wire simulations are that (i) they provide an independent 
check on the iV-body code; (ii) they isolate the effects of resonant relaxation without 
additional "noise" from non-resonant relaxation; (iii) the evolution of AL and AL can be 
traced for longer times than in the iV-body code, due to the larger timesteps that can be 
used (but see below). 

2.2. Simulation Results 

A sample pair of simulations illustrating the presence of resonant relaxation is shown in 
Figures [l] and [|, which show simulations in the Kepler and isochrone potentials respectively; 
we expect resonant relaxation to be important in the former case and not the latter. Both 
simulations used n = 4 test stars out of n + N = 64 total stars, each of mass m = 3x 10 _7 M. 
(Although the total number of stars in each cluster was small, the results they display are 
consistent with those from larger simulations; in visual terms the small- N runs show the 
effects of resonant relaxation most clearly since they can be evolved longer.) 
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Fig. 1. — Energy and angular momentum relaxation over time (in units of the mean orbital 
time, t or b) for four test stars in a Keplerian star cluster. Resonant relaxation causes AL 
(dotted line) and \AL\ (dashed line) to grow almost linearly with time, in contrast to the 
square- root growth of AE (solid line), which is consistent with a random walk. Note the 
(weak) evidence for a turnover in the AL curves at t ~ tp rec ~ 10 4 ' 7 t O rb, where resonant 
relaxation is expected to become a random walk (see eqs. |] and ||). 
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Fig. 2. — As Figure |l| but for an isochrone potential, in which there is rapid precession. 
Note that in this case AL/L ~ AE/E oc t 1 / 2 at all times; because the isochrone potential is 
spherical, however, resonant relaxation is present in \AL\ 



L4). 
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The four panels in each of Figures |T] and |2] plot the three quantities AE/E , AL/L max , 
and \AL\/L max for the four test stars, where E is the initial energy and 

L ( E ) = i GM ^ 2 ^ 1/ ^ Kepler potential; , . 

maxl ) \ GM /(2\E\) 1 / 2 (l-2b\E\/GM), isochrone potential, 1 ' 

is the angular momentum of a circular orbit with energy E. The times are given in units 
of t or b, defined as the mean radial orbital period of all stars in the cluster. The difference 
between the two potentials is clear: in the Kepler case (Fig. [I]), growth of AL systematically 
outpaces that of AE, by 1-2 orders of magnitude in the course of the integration, whereas 
in the isochrone potential (Fig. 0) AL/L max ~ AE/E at all times. In other words, if 
we define the (scalar) angular momentum relaxation time ti by AL(ti) ~ L max (£'o) and 
the energy relaxation time t E by AE{t E ) ~ E , then t L <C t E when precession is slow, 
while tL ~ t E otherwise. In both potentials, the vector angular momentum relaxation time 
tL *C t E (cf. § O); however, this rapid relaxation of L is less interesting since it reflects 



changes in the orientation rather than the shape of the orbit, which are of little consequence 
in non-rotating spherical clusters. 

Note that for these simulations, t pre c/^orb ~ M/M* « 10 4 ' 7 , so the figures mostly 
illustrate the regime t <C t prec in which AL is expected to increase linearly with time (eq. |6|). 
At times t ^> tp rec , AL should increase as the square root of the time (eq. ||). There is a 
suggestion of such a turnover in Figure [TJ; we comment further on this issue below. 

We may test the analytic analysis of § H in more detail by performing least-squares fits 
on the simulation data. We fit each of the curves in Figures and ^| — and analogous curves 
for other simulations — to the variety of power laws, sums of power laws, and broken power 
laws defined in Table |l|. Formula (a) tests whether AE{t) exhibits the t 1 ' 2 dependence 
predicted by equation (||), and formula (b) finds the best-fit value of the parameter a 



Table 1. Parametric Relaxation Curve Models' 



AE(t)/E 


AL(r)/L max 


|AX(r)|/L max 


(a) CmN 1 / 2 ^ 


(c) 7 s miV 1 /V s 


(f) 7v miV 1 /V v 


(b) amN^r 1 ' 2 


(d) r] s mN l l 2 T l l 2 + f3 s mN l ' 2 T 


(g) r/.miV^V/ 2 + (3 v mN^ 2 T 




(e) cr fel /[l + (r/Tb)^- 62 )/^ 





T = (t/t orb ); M = 1. 
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in equation (|j). Formula (c) tests whether AL(t) is dominated by the t 1//2 behavior 
characteristic of non-resonant relaxation (eq. |J) or the linear growth characteristic of 
resonant relaxation when t <C t prec (eq. []); formula (d) fits simultaneously for both resonant 
and non- resonant relaxation in AL (eqs. [| and ||); and formula (e) fits to a broken power 
law (cx t bl for small t and oc t b ' 2 for large t, the transition occurring at t/t rb ~ Tb) to look 
for the turnover from linear growth (eq. ||) to square-root growth (eq. |8|) that is expected at 
t ~ t P rec- Finally, formulae (f) and (g) are the vector analogs of (c) and (d) (cf. eqs. ^ and 

ID- 

In performing the fits, random noise in the relaxation curves of individual test stars is 
a significant concern. Note first that fitting to the raw output is undesirable because curves 
such as AE(t) cross zero many times, which in log-log space implies passage through — oo; 
what is wanted is a fit to the more steadily growing envelope, ignoring the zero crossings. 
This was done by smoothing the individual curves using an averaging window of width 
Alogt ~ 0.1 to replace each data point by the weighted average of its neighbors, larger 
data points receiving greater weight (the weighting factor was proportional to the square 
of the data point). This window was large enough to remove the sharp "dropouts" due 
to zero crossings, but too small to affect the global shape of the envelope (decreasing the 
window size to Alogt ~ 0.03 changed power-law slopes by ^ 1%, but noticeably increased 
the random noise). As these smoothed curves still contain (real!) random variations, 
more robust results were obtained by averaging the curves for all test stars in a particular 
simulation, scaled temporally by their respective orbital periods, to create composite average 
relaxation curves for the cluster, model parameters being calculated from these composites. 
Conceptually one can think of these curves as belonging to a fictitious "average" star with 
orbital period t or b; we will use this term to label such fits. An alternative procedure is to 
derive best-fit parameters for individual stars and then average these results; this method 
is less robust since it attempts to remove the random noise after performing the fits, which 
can be ineffective if this noise causes the fits to be poor. Therefore, all of our quantitative 
results are based solely on fits to the average star as defined above. 

Tables H and |3| show derived parameter values for the relaxation curves shown in 
Figures [l] and ||. Fits to the individual test stars, the averages of the fit parameters (labeled 
mean), and fits to the average star (labeled AS) are given. Comparison of Tables |2| and 
|3] shows the effect of resonant relaxation quite clearly. In Table |2] (Kepler model) the 
exponent 5 S for AL/L mSuX (formula [c]) is near unity, corresponding to the linear growth 
expected for resonant relaxation; while the exponent v for AE/E (formula [a]) is near 0.5, 
corresponding to a random walk dominated by non-resonant relaxation (recall that there is 
no resonant energy relaxation). On the other hand, in Table |3] (isochrone model) there is 
no systematic difference between the exponents S s and v. 
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Table 2. Best-fit Parameters for Figure [I] (Kepler Potential) 



Formula Number and Parameters 





(b) 


(a) 


(c) 


(f) 


(d) 


(g) 




(e) 


Star 


a 


C 


v 


7s, h 


7v, K 


Ps 


Vv, Pv 


log 


r b , b 2 


1 


1.73 


3.8 


0.39 


2.0, 0.92 


2.4, 0.96 


2.8, 0.92 


3.2, 1.51 


4.7 


-0.17 


2 


5.40 


9.9 


0.40 


0.9, 0.84 


3.4, 0.94 


1.2, 0.22 


2.5, 2.03 


3.8 


1.05 


3 


2.31 


2.9 


0.46 


2.2, 1.00 


2.5, 1.00 


0.0, 2.50 


0.0, 2.76 


4.5 


-0.50 


4 


6.36 


5.2 


0.53 


2.4, 0.88 


5.6, 0.87 


0.6, 1.17 


4.0, 2.20 


3.8 


-0.10 


mean 


3.26 


4.2 


0.45 


2.0, 0.95 


3.2, 0.95 


1.9, 1.28 


2.2, 2.10 


4.4 


0.07 


AS 


3.51 


5.2 


0.44 


1.8, 0.90 


3.5, 0.93 


1.6, 0.85 


2.4, 2.04 


4.5 


0.18 



Table 3. Best-fit Parameters for Figure ^| (Isochrone Potential) 



Formula Number and Parameters 





(b) 


( 


») 


(c 





(f) 


(d) 


(g) 


(e) 


Star 


a 


C 


V 


7s, 


6 S 


7v, <5 V 


Vb, Pb 


Vv, Pv 


logr b , b 2 


1 


2.36 


1.7, 


0.56 


2.3, 


0.52 


2.3, 0.85 


2.5, 0.00 


3.7, 0.55 


0.1, 0.52 


2 


5.06 


3.0, 


0.59 


2.0, 


0.61 


1.8, 0.91 


3.1, 0.01 


2.1, 0.75 


-0.1, 0.61 


3 


2.87 


8.5, 


0.34 


7.0, 


0.34 


4.4, 0.89 


3.4, 0.00 


6.5, 1.47 


0.6, 0.33 


4 


1.65 


2.7, 


0.43 


3.1, 


0.44 


1.5, 0.95 


2.3, 0.00 


2.6, 0.82 


0.5, 0.44 


mean 


2.80 


3.5, 


0.48 


2.8, 


0.50 


3.0, 0.89 


2.6, 0.00 


4.6, 0.99 


0.2, 0.49 


AS 


2.84 


3.2, 


0.48 


3.1, 


0.48 


2.5, 0.89 


2.9, 0.00 


3.6, 0.84 


0.5, 0.46 
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The difference between the Kepler and isochrone potentials is also reffected in the fits 
to formula (d), which decomposes the relaxation into terms oc t 1 ^ 2 (non-resonant) and oc t 
(resonant). The strength of the resonant term, which is proportional to /3 S , is comparable 
to the strength of the non-resonant term in the Kepler potential, but is essentially zero in 
the isochrone model. 

The presence of resonant relaxation of the vector angular momentum in spherical 
potentials is illustrated by the exponent 5 V for AL/L max , which is near unity in both 
Tables and |(eq. [15]). 

A glance at the derived values for log and 62 in Table illustrates the advantage 
of using the average star to compute parameters instead of averaging the fits to 
individual stars; the individual fits are so noisy that averaging after the fact is nearly 
meaningless. Note that log Tb for the average star in Table || is of order the expected value, 

log(VecAorb) ~ l0g(M/M*) ~ 4.7. 

We originally hoped that the iV-wire code could be run for many more orbital times 
than the corresponding iV-body calculation, since the required integration step is much 
longer. This hope was frustrated by close encounters of the wires, which required each wire 
to be divided into hundreds or even thousands of segments for accurate calculation of the 
torque (the calculations for each segment taking ~ 10 2 -10 3 times longer than the simple 
force calculations of the iV-body code). These encounters dominated the computation 
time, so that the iV-wire code was significantly slower than the iV-body code, particularly 
for larger N. Nevertheless, the iV-wire simulations that were performed showed the same 
qualitative features as their iV-body counterparts, including clear evidence of resonant 
relaxation of similar strength, and the suggestion of a turnover near t ~ t pTec . The iV-wire 
code could be sped up in several ways: by softening the potential from the wires, by using 
a finer mesh close to the crossing point, or by analytic evaluation of the torque when 
the wires are close to crossing. However, as there appeared to be no new features in the 
relaxation curves, a substantial effort to improve the speed was not deemed worthwhile, 
and quantitative results from the wire simulations will not be given. 

Table |4] lists the best-fit model parameters averaged over « 30 iV-body simulations 
spanning several orders of magnitude in m (from 3 x 10 -8 to 3 x 10~ 5 ) and iV (from 64 to 
8192). The contribution of each simulation to the average was weighted by the standard 
error in the fit, which in turn was determined by choosing the error for each data point in 
the simulation so that x 2 — 1 P er degree of freedom. 

The results nicely confirm our theoretical expectations: (i) In both the Kepler and 
isochrone potentials, the evolution of AE is consistent with a random walk, that is, the 
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Table 4. Global Average Best-fit Parameters 3, 



Formula 


Parameters 


Keplei 


- Potential 




Isochrone Potential 


(b) 


a 


3.09(0.11) 




[0.25] 


2.76(0.11) 


[0.23] 


(a) 




3.24(0.17), 


0.49(0.01) 


[0.23] 


2.93(0.18), 0.49(0.01) 


[0.21] 


(c) 


7s, S 8 


1.35(0.08), 


0.85(0.02) 


[0.26] 


2.87(0.12), 0.49(0.01) 


[0.21] 


(f) 


7v, 5v 


2.78(0.18), 


0.88(0.01) 


[0.17] 


2.73(0.15), 0.84(0.02) 


[0.33] 


(d) 


?7s, /5s 


1.37(0.11), 


0.53(0.06) 


[0.22] 


2.76(0.09), 0.00(0.01) 


[0.24] 


(g) 




2.08(0.26), 


1.79(0.12) 


[0.11] 


3.50(0.10), 0.71(0.07) 


[0.13] 


(e) 


logr b , 6 2 


0.02(0.27) b , 0.49(0.13) 


[0.04] 


0.73(0.14), 0.49(0.01) 


[0.04] 



a The numbers in parentheses are 1-a errors on the best-fit parameters, which were derived 
by assuming equal error (in dex) for each simulation data point and choosing this error so 
that x 2 = 1 per degree of freedom (the required errors in the data points in dex are given in 
square brackets). 

b This number is log(rb/r prcc ), not logrb. 
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exponent v = 0.5 to within the errors; moreover, the rate of energy relaxation is nearly 
identical in the two models (the values for a differ by only 10%), which confirms that they 
are close analogs except for the resonance in the Kepler case, (ii) In the isochrone potential, 
the evolution of AL is consistent with a random walk, that is, 5 S = 0.5 to within the errors; 
while in the Kepler case AL grows much more rapidly (the derived value 5 S = 0.85 ± 0.02 is 
less than the expected value of 1, presumably because small-scale fluctuations in the torque 
add a random- walk component — this issue is discussed further below). There is also no 
evidence for a linear component in AL in the isochrone model (/3 S = 0). (iii) The growth 
of AL is more rapid than a random walk in both potentials; once again, the derived values 
5 V = 0.88 ± 0.01, 0.84 ± 0.02 are lower than unity because of small-scale fluctuations in 
the torque, (iv) The value of r b is close to the expected value r prec in the Kepler potential 
(log(r b /r prec ) = 0.02 ± 0.27); moreover, as predicted by equation (^), the evolution of AL 
for t tp r ec is consistent with a random walk (which requires b<i = 0.5; the observed value 
is 0.49 ±0.13). 

Note that the diffusion in angular momentum, as measured by r] s and r] v , is stronger 
in the isochrone than the Kepler potential (perhaps because the encounter velocities are 
smaller in the isochrone core), while the resonant relaxation in L, as measured by (3 V , 
is stronger in the Kepler potential because of the contribution to the torque from the 
eccentricity of the orbits. 

We now comment on whether the linear /random- walk model (formulae [d] and [g]) 
provides a better fit to AL and AL than a simple power-law (formulae [c] and [f]). Table ^ 
shows that in all three cases of interest (Kepler AL and AL, and isochrone AL) formula 
(d) fits better than (c) and (g) fits better than (f), the required errors in the data points 
to achieve a given x 2 being smaller. To determine whether this difference is significant, 
paired sample Student's t-tests were performed to compute the confidence level at which 
the hypothesis of the corresponding x 2 values having the same mean could be excluded. 
For both potential models, the resulting confidence for AL was over 99.9% — the power-law 
model is clearly inferior to the linear /random- walk model. In the case of AL (Kepler 
potential), applying the t-test to the full simulation set gave a confidence of only « 85% 
that the linear /random- walk model provides a better fit; this weaker result is not surprising 
since the relaxation curves are noisier. By performing the test using only the 6 simulations 
with at least 32 test stars, we obtained a stronger confidence level of ~ 95%. 

We also checked the scaling of the relaxation rate with m and N, by fitting the 
simulations to a formula of the form AE/Eq = am a iVV/ 2 ; this gave a = 3.6 ± 0.9, 
a = 1.02 ±0.01, and b = 0.51±0.03. Thus the exponents are consistent with their predicted 
values a = 1, b = ~ to within a few percent. Similar results were obtained from fits to AL 
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and AL. 

We note that softening is expected to affect the strength of the random walk terms 
in the simulations by changing the value of the "In A" factor which these terms implicitly 
contain. In general, one expects A ~ Rj ma.x(s,GM/v 2 ), where s is the softening length 
and R the system size. Since s 3> GM/v 2 in all of our simulations, we expect that 
a 2 oc t^ 1 oc hi[K (R/s)] for some factor K ~ 1. To test this scaling, we ran a series of 
Kepler simulations in which the softening length was reduced by a factor of 100 (from 
s ~ 10~ 2 -R to s ~ 10 _4 -R), which should increase a by a factor of about 1.4. It was found 
that a = 5.5 ± 0.2 for these runs, a factor of 1.8 larger than in Table |], which is in line with 
expectations given the uncertainties. Similar increases were seen in other parameters, such 
as rj s which increased by a factor of 1.6. This observation also supports our claim that the 
7] parameter really is measuring a random-walk component in AL. 



3. Implications for Galactic Nuclei 

3.1. Black Hole Fueling Rates 

Tidally disrupted stars can be an important fuel source for massive black holes in 
stellar systems, disruption occurring if the stars pass closer to the black hole than their 
tidal radius, r t . The pericenter of an orbit with energy E ^> —GM/r t is within the tidal 
radius if its angular momentum is less than L m j n c±. (2GMr t ) 1//2 ; the region in phase space 
L < L m in is known as the "loss cone" because stars in this region will be disrupted in 
less than one orbital period, unless their orbits are perturbed out of the loss cone before 
reaching pericenter. Angular momentum relaxation provides a steady supply of stars to 
the loss cone. This mass supply rate is interesting in the context of active galactic nuclei 
(AGNs), for which tidal disruption of stars is a possible fueling mechanism. Previous 
studies (e.g., Duncan & Shapiro 1983; David et al. 1987a,b; Murphy et al. 1991; Polnarev 
& Rees 1994) concluded that (non-resonant) relaxation is too slow to sustain typical AGN 
luminosities, unless mass loss is dominated by other mechanisms (such as physical collisions 
between stars) or additional, large-scale torques are present (as from a nuclear bar, binary 
black hole, or triaxial galactic potential). In this section we investigate whether resonant 
relaxation can enhance the supply of stars to the loss cone and hence the fueling rate of 
AGNs from disrupted stars. 

We shall parameterize the loss rate of stars by the dimensionless number A(L), which 
is the fraction of stars at a given energy that are disrupted each orbital period. Let AL or b 
be the root-mean-square change in scalar angular momentum in one orbital period. For 
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non-resonant relaxation, equation (W) implies 



^ ~ vJ^r, (38) 
while for resonant relaxation, equation (Q) implies 

the non-resonant contribution to AL or b is stronger by of order the square root of the 
Coulomb logarithm. Lightman & Shapiro (1977) distinguish two cases: the "pinhole" limit 
in which AL or b 3> L min , and the "diffusion" limit in which AL or b <C L min . First consider the 
pinhole limit. Here leakage of stars into the loss cone is too small to affect the distribution 
function significantly, so that the distribution function is approximately uniform on the 
energy surface in phase space, even near the loss cone. Hence the fraction of stars at a given 
energy that is lost per orbital period is simply the fractional area in the energy surface 
occupied by the loss cone, 

X{E) ^1TJE) ( AL orb^^nin). (40) 

11 1 ctX \ / 

In this limit the fueling rate is independent of the strength of the relaxation, so long as 
AL or b ^ Lm in - Therefore the loss rate in this case is unaffected by resonant relaxation, since 
AL or b is dominated by non-resonant relaxation. The pinhole limit requires 



& ■ (41) 



L max M 



When is replaced by "=," equation pi] ) defines the so-called critical radius, r crit ; in 
general the pinhole limit applies outside the critical radius and the diffusion limit applies at 
smaller radii. 

Next consider the diffusion limit. In the absence of resonant relaxation the loss rate 
can be determined by solving the Fokker-Planck equation (e.g., Lightman & Shapiro 1977; 
Cohn & Kulsrud 1978), and is found to be 

^(E) * T 2 rk\ (r IT \ ( AL - b £ W ( 4 2) 

Now assume resonant relaxation is present. Over timescales <C t prec the resonant torque is 
approximately constant, so the angular momentum drifts at a nearly uniform rate. Since 
^precAorb ~ M / M± 3> 1 whenever resonant relaxation is important, we may assume that 
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(^prec/^orb) AL™ b ^ Anhu at least if we are not too far inside the critical radius. Then the 
evolution of the angular momentum on the scale of the loss cone resembles a steady drift 
rather than diffusion, so a typical area ~ L min AL^ is swept into the loss cone per orbital 
period, and the loss rate is 

T . A T rcs 

A res (E) ~ ™ orb ( AL orb £ L min ) . (43) 

III cLX \ / 

The ratio of the resonant loss rate to the non-resonant loss rate in the diffusion limit is 

\ res T A j res 

^ ~ TXT^f 1 M WW ( AL orb £ L min ) . (44) 
Since AL ^ b and AL™ b differ only by a logarithmic factor (cf. eqs. |38| and f39fl), the loss 



rate in the diffusion limit is much larger under resonant relaxation than under non-resonant 



relaxation. (Note that formula [53J neglects the relativistic effects discussed in § |3.2| , which 
reduce the resonant loss rate in galactic nuclei.) However, this enhanced loss rate does 
not strongly affect the total rate of fueling of the black hole by bound stars. The reason 
is that the overall flux of stars into the cusp — which in a steady state equals the fueling 
rate from bound stars — is determined by the bottleneck at the critical radius (i.e., by the 
slow diffusion of orbital energies, the rate of which is unaffected by resonant relaxation); 
the enhanced loss rate inside r crit will reduce the density of stars at r < r crit and the radial 
profile of the cusp of bound stars, but it will not strongly affect the overall disruption rate. 
We conclude that resonant relaxation does not significantly enhance the viability of tidal 
disruption as a possible AGN fueling mechanism. 



3.2. Relativistic Effects 

We have assumed so far that stellar orbits are Keplerian whenever the black-hole mass 
is much larger than the mass in stars. At small radii, however, relativistic effects cause orbits 
to precess rapidly, which quenches the resonant relaxation. In the present case the most 
important relativistic effects are the precession of periapsis and Lense-Thirring precession 
of the orbital plane. The former effect (observed in the solar system as a component 
of Mercury's perihelion precession, for example) is present around both Schwarzschild 
(non-rotating) and Kerr (rotating) black holes; in both cases, the change in the argument 
of periapsis per orbit in the weak-field limit is given to leading order by 



(45) 
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where M is the hole mass, a and e are the orbit's semi-major axis and eccentricity, and 
the precession is in the direction of motion (e.g., Darwin 1961). Around a Kerr black 
hole — the expected case for the nuclei of AGNs, since disk accretion will spin up an initially 
non-rotating black hole (Bardeen 1970, Thorne 1974) — orbits also undergo Lense-Thirring 
precession, which in the weak-field limit appears as a precession of the orbital angular 
momentum vector around the black hole's angular momentum vector. However, this effect 
merely adds another (uninteresting) secular growth term to AL, and hence is not important 
for the present analysis. 

Let us define a relativistic precession timescale, t GR c = [tt/Alo] t or b, and a non- 
relativistic precession timescale, t^ c ~ (M/M*)t orb (eq. [I]), where t or b is the orbital 
period. Relativistic precession can be important only when t GR c ^ t^ c . m fact, since 
the relativistic precession is prograde (Au; GR > 0) while the non-relativistic precession is 
generally retrograde (Auj Kcp < 0), relativistic effects slightly enhance resonant relaxation 
so long as t GR c ^ t^ c . Even when t GR c < t^ c , resonant angular momentum relaxation 
remains important and is described by equation (^) so long as i GR c ^> t or b, although the 
relaxation time is increased by a factor tp^ c /t GR c . 

The angular momentum of a Kepler orbit is given by L 2 = GMa(l — e 2 ), so the 
relativistic precession time may be written 

«■ " i (m) ** (46) 

Thus resonant relaxation is completely quenched by relativistic effects when L <^ L GR , where 
L GR = GM/c. A more precise statement is that there is a region of angular- momentum 
space in which resonant relaxation is completely quenched only if 

L GR £ max(AL orb , L min ); (47) 

if L GH < L m i n stars are disrupted before the relativistic precession time becomes comparable 
to the orbital time, and if L GR < AL or b a star can relax across the relativistic region 
in less than one orbit, thereby avoiding the relativistic precession (nearly all of which 
occurs at periapsis). If the inequality ( ^7|) is satisfied, the angular momentum diffuses 
through resonant relaxation when L ^ L GR and by non-resonant relaxation when 
L min + AL or b ^ L ^ L GR ; stars will pile up in the latter region of phase space since 
non-resonant relaxation is slower. Even outside this region, relativistic precession can 
dominate the precession rate and hence reduce the rate of resonant relaxation. 

The inequality fl47j) requires either that the loss cone is empty (AL orb ^ L min ) and 
that L GR /L min ^ 1, or that L GR /AL or b ^ 1 if the loss cone is full; since resonant relaxation 
does not affect disruption rates if the loss cone is full, the latter case is uninteresting. A 
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numerical value for the former condition is easily derived; noting that L min ~ (2GMr t ) 1 / 2 , 
where r t ~ 2r*(M/m) 1 / 3 is the tidal radius of a star of mass m and radius r*, assuming 
r* oc m 2 / 3 (reasonable for both low and high mass main-sequence stars), and normalizing to 
solar values, we obtain 

L GR __/ r M \ 1/3 /M \ 1/6 



Resonant relaxation only dominates non-resonant relaxation when t prcc ^ 10t or b (cf. 



eq. |52D — which here implies L ^ 8L GR . Therefore we may assume that relativistic precession 
will substantially affect the loss rate when the loss cone is empty and L min ^ 4L GR , which 
implies 

/ m \ 1/2 

m ^ Ax1 ° 7 {m^) Mq - (49) 

The mass estimate (|49"D is quite uncertain, but lies squarely inside an interesting range: the 
estimated masses of black holes in galactic nuclei range from ~ 10 6 Mq to over 10 9 Mq, 
and the Eddington luminosity corresponding to fl49|) is LecM ~ 5 x 10 45 erg s _1 — a typical 
AGN luminosity. 

To summarize, in galactic nuclei relativistic precession substantially reduces — and can 
even completely quench — resonant relaxation close to the loss cone. A more sophisticated 



treatment than the one leading to equation fl43 ) would be required to determine the 



disruption rate and density profile when the loss cone is empty. Such an analysis is beyond 
the goals of this investigation. 



3.3. Relaxation Time Estimates 

The non-resonant relaxation time in a stellar system may be written (e.g., Binney & 
Tremaine 1987, eq. 8-71) 

C = 0-3 ^ 2 r , (50) 
G z mp in A 

where a is the one-dimensional velocity dispersion and p is the stellar density. Specializing to 
the case of stars in the potential of a black hole of mass M, we may write a ~ (GM/3?") 1 / 2 , 
p ~ M*/(|-7rr 3 ), and A ~ A; thus 

jm-3/2 3/2 

Ci * 0.3- 



G l l 2 m 2 NhiN 



4 x 10 10 yr / M x / ™ \ x / 2 /i /i^ x / ^ \ 3 / 2 



In A \MJ \10 S M Q ) V m / \ 1 pc J ' ( ' J ''' ) 



-28- 



The resonant relaxation time (cf. §§ 1.2- |1.3|) is 



n 2 M M 

f res ^ V^^* ,nr ^ j ^* ,nr (coN 

* rel " # M * rel - 7 M * re1 ' (52j 

where /3 S and r? s are taken from Table [|. 

Table [5] applies these estimates to several nearby galaxies that are believed to harbor 
massive black holes (Kormendy & Richstone 1995). The table lists the resonant and 
non- resonant relaxation times t r r lf and t™ h at O'.'l from the center for external galaxies 
(roughly Hubble Space Telescope resolution) and at l'.'O for the Galaxy. 

We see that resonant relaxation generally does not influence the structure of galaxies 
(^-rei ^ 10 10 y r ) & t available resolutions, except in a few Local Group members such as M32 
and the Galaxy. M31 contains a region in which tUf < 10 10 yr < tf^ that is almost (!) 
resolvable by HST. At smaller distances from the black hole the resonant relaxation time 
is shorter; crude extrapolations indicate that t^i ^ 10 9 J T a1, radii less than 0'.'05 in M32, 
at O'.'Ol in M31, and at 0"002 in NGC 3115. The first two of these could be resolved with 
proposed space-based interferometers. 



4. Discussion 

Resonant relaxation is the dominant source of angular momentum relaxation for stellar 
systems in near-Keplerian potentials, and thus plays an important role in determining the 
structure of stellar cusps around black holes in galactic nuclei or globular clusters. Resonant 
relaxation enhances the angular momentum relaxation rate by roughly the ratio of the mass 
of the black hole to the mass in stars but does not affect the energy relaxation rate (more 
precisely, the combination of actions J± — J2 is conserved, where J\ is the radial action and 
J2 is the angular momentum). 

Resonant relaxation is also present in the harmonic potentials that characterize 
constant-density cores, and may enhance the rate of angular momentum relaxation in the 
cores of globular clusters. In const ant- density cores, resonant relaxation preserves the 
combination of actions J\ — 2J2. This form of resonance is not likely to be important for 
elliptical galaxies, which do not generally have constant-density cores (e.g., Gebhardt et al. 
1996). 

One might speculate that generic potentials contain (high-order) resonances that 
are strong enough to support resonant relaxation. In this case the angular momentum 
relaxation time would be much shorter than the energy relaxation time throughout most 
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Table 5. Characteristics of Nearby Massive Black Hole Candidates 



Galaxy 


D 


M 


6(fi = 0.1) 


/ u(0'.'l) a 


^(o"i) a 


Ci(0"l) a 


a 




(Mpc) 


(M Q ) 


(") 




(yr) 


(y r ) 


CO 


Milky Way 


0.0085 


2 x 10 6 


6 


0.001 


2 x 10 s 


3 x 10 10 


12 


M32 


0.7 


2 x 10 6 


0.08 


0.2 


4 x 10 9 


3 x 10 9 


0.1 


M31 


0.7 


3 x 10 7 


0.5 


0.002 


2 x 10 10 


1 x 10 12 


0.07 


NGC 3377 


9.9 


8 x 10 7 


0.05 


0.3 


1 x 10 12 


3 x 10 11 


0.005 


NGC 3115 


8.4 


2 x 10 9 


0.5 


0.02 


2 x 10 12 


2 x 10 13 


0.002 


M87 


15.3 


3 x 10 9 


1.3 


0.001 


8 x 10 12 


1 x 10 15 


0.001 



For the Galactic Center only, these values are given at l'.'O, not O'/l. 



Note. — Estimated distance D and black-hole mass M are from Kormendy & Richstone 
(1995). The quantity fi(r) is the fraction of the total mass inside radius r that resides in stars; 
/i < 1 implies that the potential is nearly Keplerian and when [i <^ 0.1 resonant relaxation 
dominates non- resonant relaxation (eq. |52|). The resonant and non-resonant relaxation times 
t T T lf and t"ei are given by equations (^1|) and fl5"2|) ; 9 Ie8 is the apparent angular distance from 
the black hole at which the resonant relaxation time is 10 10 yr. Stellar density estimates 
for external galaxies are based on Hubble Space Telescope photometry (Faber et al. 1996) 
and assumed mass-to-light ratios; densities for the Milky Way are based on near-infrared 
photometry (Kent 1992) and an assumed core radius of 0.15 pc (Eckart et al. 1993). 
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of a galaxy. We suspect that this speculation is not correct, since our iV-body simulations 
in the isochrone potential yield very similar relaxation rates for the energy and angular 
momentum. Nevertheless, the presence of resonant relaxation is a reminder that our 
understanding of relaxation in stellar systems is crude, and has not been numerically 
verified under conditions (N ~ 10 11 ) found in real galaxies. 

Resonant friction leads to growth or decay of the eccentricity of massive objects 
orbiting in near-Kepler potentials, depending on whether the star orbits are predominantly 
radial or tangential. Resonant friction can strongly influence the orbital evolution of a 
binary black hole (at least if the mass ratio of the binary is sufficiently far from unity). 
In radially biased star clusters the eccentricity of the binary will grow, at a rate faster 
than the decay of the orbital energy, at least if the friction is dominated by cluster stars 
rather than unbound stars. The binary eccentricity will grow until the resonant friction 
is quenched by relativistic precession, at which point gravitational radiation may erode 
the energy of the binary faster than non-resonant dynamical friction. The details of this 
evolution are relevant to the merger rate of black holes, the gravitational-wave background, 
the prevalence of binary black holes in AGNs, and the viability of massive black holes as 
dark matter candidates (see Quinlan 1996 for references). 

Resonant friction can also erode the inclination of a massive object in a rotating stellar 
system. This process may be relevant to a star cluster in which there is a massive accretion 
disk (Ostriker 1983; Syer et al. 1991); resonant friction could accelerate the evolution of 
massive stars into low-inclination orbits embedded in the accretion disk. 

The analytical treatment of resonant relaxation that we have offered in § |1] is only 
approximate. Accurate expressions for the resonant and non-resonant relaxation rates in 
a given star cluster could be derived by expanding the potential from a stellar orbit in 
action-angle variables. So far this procedure has only been carried out for the dynamical 
friction component (Weinberg 1986). The relative simplicity of the diffusion coefficients that 
describe non-resonant relaxation (e.g., Binney & Tremaine 1987) is illusory in near-Keplerian 
and other near-resonant potentials — except to describe energy relaxation — since resonant 
relaxation is stronger, and depends more sensitively on the structure of the stellar system. 
For order-of-magnitude estimates we have used the formulae in § [I], with the dimensionless 
coefficients given in Table [§. 



The estimates of tidal disruption rates in § |3.1| suffer from the absence of a consistent 
treatment of relativistic precession, which detunes the Kepler resonance near the loss cone 
in galactic nuclei. However, resonant relaxation is unlikely to increase substantially the 
tidal disruption rate, which is mostly determined by the location of the critical radius 
r crit , set by the angular momentum changes in a single orbital period. For similar reasons, 
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resonant relaxation will not greatly affect the disruption rate resulting when the wandering 
("Brownian motion") of the black hole from the center of the nucleus (Quinlan 1995) 
is taken into account. Thus tidal disruption appears to remain incapable of powering 
typical AGNs; however, since the relativistic detuning is largely ineffective at hole masses 
<^ 10 7 Mq (eq. |48"D , resonant effects may offer modest improvements in the feasibility of 
disruption-dominated mass loss in Seyferts and other nuclei containing low-level AGN 
activity, for which the energy requirements are less severe. Similarly, resonant relaxation 
may modestly enhance the rate of flares from tidally disrupted stars in nearby galaxies with 
central black holes (Rees 1988). 

The effectiveness of relativistic precession in disabling resonant relaxation illustrates 
that general relativity can have dramatic physical consequences even where the motion is 
predominantly Newtonian; in particular, the shape of the density cusp inside r cr i t can be 
strongly dependent on relativistic precession. Thus resonant relaxation might one day be 
used to show that the massive dark objects observed in galactic nuclei are indeed black 
holes (or at least behave as such on a scale of ~ 10 2 Schwarzschild radii) — a conclusion 
which today must be reached by indirect (albeit compelling) arguments. 

The discussion in §3.1 also illustrates that the Fokker-Planck equation used to describe 
non-resonant relaxation (e.g. Binney and Tremaine 1987) is not always adequate to describe 
resonant relaxation. The Fokker-Planck equation assumes that the fluctuating forces at 
different times and locations are uncorrelated, i.e., that the correlation function of §1.1 
has the form Cy(ri, V2, t) = Kij(ri)5(r2 — ri)5(r). This is a reasonable approximation 
for non-resonant relaxation, which is dominated by close encounters (cf. §1.1). In 
contrast, the resonant forces are correlated over large spatial scales and over times ~ t prec . 
The inadequacy of the Fokker-Planck approximation (or the master equation, or the 
approximation that relaxation is a Markov process), is particularly acute in the diffusion 
limit (eq. f£|), when the size of the loss cone in angular-momentum much 
greater than the change in angular momentum per orbit AL or b but much less than the 
change in angular momentum per precession time. 

There is an appealing analogy between relaxation of stars in angular-momentum space 
and models of stellar structure. Non-resonant relaxation is a random walk in L-space, as for 
the motion of ions in the radiative zone of a star. Resonant relaxation implies a large-scale 
drift in L superimposed upon a small-scale random walk, analogous to ionic motion in a 
convective stellar envelope. The quenching of resonant relaxation by relativistic precession 
can produce a random-walk dominated "core" in L-space surrounded by a drift-dominated 
envelope — conceptually similar to the radiative core/convective envelope structure found in 
solar-type stars. There are fundamental differences: L-space has no analog to gravity, but 
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there is a net flux of stars towards small L due to removal of stars by tidal disruption. 

The structure of a relaxed star cluster surrounding a black hole has been examined by 
several authors (Peebles 1972; Bahcall & Wolf 1976, 1977; Lightman & Shapiro 1977; Cohn 
& Kulsrud 1978). These analyses do not take resonant relaxation into account and therefore 
some of their conclusions may be suspect: we expect that including resonant relaxation will 
not strongly affect the structure of the star cluster outside the critical radius r crit (§ |3.1| ) or 
the total flux of stars into the loss cone, but may substantially reduce the density of stars 
inside r cr i t . This classic problem should be re- investigated. 

Resonant relaxation implies that there may be regions near the centers of elliptical 
galaxies (typically <^ 1 pc in radial extent; cf. Tabled) that are relaxed in angular momentum 
but not energy. If non-rotating, such regions will have isotropic distribution functions; if 
rotating, the mean rotation speed will depend on the stellar mass. Unfortunately, these 
regions are not accessible at Hubble Space Telescope resolutions in most nearby galaxies. 
A more fundamental problem regarding the possible observational detection of resonant 
relaxation is that the isotropy it produces will be undetectable unless the initial distribution 
function is significantly anisotropic; there is currently no evidence for this in observed 
nuclear star clusters. 

We thank Jihad Touma for supplying the program to compute the average torque 
between two Keplerian orbits. This research was supported by NSERC, and by a Jeffrey L. 
Bishop Fellowship to K. R. 
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